!! Module to be used if particle present
module mod_part_interp
use omplib
use technical
use general
implicit none
save
!
double precision,allocatable,dimension(:,:)::Xp,Yp,Vxp,Vyp,Uxp,Uyp,Axp,Ayp
double precision,allocatable,dimension(:,:)::Lamp,dxuxp,dxuyp,dyuxp
integer:: Nparticle,nst,nprun
double precision :: deltax,deltay
double precision,allocatable,dimension(:)::tau,one_over_tau
logical :: lcalc_Gup
logical :: write_fullparticle
!
private
public :: nprun,write_fullparticle
public :: allocate_particles,init_particle,evolve_particles,get_particle_rhs
public :: lambda_at_particle_position
public :: calc_particle_diagnostic,get_particle_coordinate
public :: read_particles_param
public :: read_input_particle,open_output_particle,write_output_particle,write_lambda_particle, &
          write_trajectory_particle,close_output_particle
namelist /particle_pars/ &
  Nparticle,nprun,nst,write_fullparticle
contains
!! ---------------SUBROUINES ---------------
subroutine read_particles_param(unit,ierr)
  integer, intent(in) :: unit
  integer, intent(inout), optional :: ierr
!----------
  if (present(ierr)) then
    read(unit,NML=particle_pars,IOSTAT=ierr)
  else
    read(unit,NML=particle_pars)
  endif
  rewind(unit)
endsubroutine read_particles_param
!!-------------------------------------------
subroutine  allocate_particles(Np,nstokes,lcalc_gradientp,allocated_memory)
  integer :: Np,nstokes
  logical :: lcalc_gradientp
  double precision :: allocated_memory
  if(Np.eq.0) then
    Nparticle=0
    nst=0
    allocated_memory=allocated_memory+0.
    lcalc_Gup=.false.
  else
   call alloc_particles(Np,nstokes,lcalc_gradientp,allocated_memory)
  endif
endsubroutine  allocate_particles
!! ----------------------------------------------
subroutine alloc_particles(Np,nstokes,lcalc_gradientp,allocated_memory)
  integer :: Np,nstokes
  integer :: ist
  logical :: lcalc_gradientp
  double precision :: allocated_memory
  Nparticle=Np
  nst=nstokes
  lcalc_Gup=lcalc_gradientp
  if (nst.gt.0) then
    allocate(tau(nst-1),one_over_tau(nst-1))
    open(unit=11,file='tau_st.in',status='old')
    do ist =1,nst-1  
      read(11,*)tau(ist)
      if (tau(ist) .gt. 1e-16) then
        one_over_tau(ist)=1./tau(ist)
      else
        write(*,*) 'WARNING: stopping time too small'
      endif
    enddo
    close(11)
    write(*,*)'The stopping time of the particles are:'
    write(*,*) tau
  endif
  allocate(Xp(Np,nstokes),Yp(Np,nstokes))
  allocated_memory=allocated_memory+2.*8.*dfloat(Np*nstokes)
  Xp=0.;Yp=0.;
  allocate(Vxp(Np,nstokes),Vyp(Np,nstokes))
  allocated_memory=allocated_memory+2.*8.*dfloat(Np*nstokes)
  Vxp=0.;Vyp=0.
  allocate(Axp(Np,nstokes),Ayp(Np,nstokes))
  allocated_memory=allocated_memory+2.*8.*dfloat(Np*nstokes)
  Axp=0.;Ayp=0.
  allocate(Uxp(Np,nstokes),Uyp(Np,nstokes))
  allocated_memory=allocated_memory+2.*8.*dfloat(Np*nstokes)
  Uxp=0.;Uyp=0.
  if (lcalc_gradientp) then
    allocate(Lamp(Np,nstokes))
    allocated_memory=allocated_memory+8.*dfloat(Np*nstokes)
    Lamp=0.
    allocate(dxuxp(Np,nstokes),dxuyp(Np,nstokes),dyuxp(Np,nstokes))
    allocated_memory=allocated_memory+3.*8.*dfloat(Np*nstokes)
    dxuxp=0.; dxuyp=0.; dxuxp=0.
  endif
!
endsubroutine alloc_particles
!%--------------------------------------
subroutine init_particle(n1,n2,dx,dy,ux,uy)
  implicit none
  integer :: n1,n2
  double precision :: dx,dy
  double precision, dimension(n1,n2) :: ux,uy
  integer::i1,st,intx,inty
!
  deltax=dx;deltay=dy
!! Initialize particle position
  call random_seed()
  call random_number(Xp(:,1)); call random_number(Yp(:,1))
!!
  do i1 = 1,Nparticle
    intx = floor(Xp(i1,1)*dfloat(n1))+1
    inty = floor(Yp(i1,1)*dfloat(n2))+1
    Xp(i1,:) = intx*dx; Yp(i1,:) = inty*dy;
    Vxp(i1,:)=ux(intx,inty)
    Vyp(i1,:)=uy(intx,inty)
  enddo
!!----------------
end subroutine init_particle
!%---------------------------------
subroutine evolve_particles_euler(dt)
  double precision :: dt
  integer :: ist
!!
!
!The first set of particles is always Lagrangian
!
!$OMP WORKSHARE
  Xp=Xp+Vxp*dt
  Yp=Yp+Vyp*dt
!$OMP END WORKSHARE
  do ist=2,nst
!$OMP WORKSHARE
    Vxp(:,ist)=Vxp(:,ist)+Axp(:,ist)*dt
    Vyp(:,ist)=Vyp(:,ist)+Ayp(:,ist)*dt
!$OMP END WORKSHARE
  enddo
!
endsubroutine evolve_particles_euler
!!----------------------------------------------
subroutine get_particle_rhs(ux,uy,n1,n2)
  integer :: n1,n2
  double precision,dimension(n1,n2) :: ux,uy
  integer :: ist
!
  call get_velocity_at_particle_position(ux,uy,n1,n2)
!
!For Lagrangian particles, particle velocity is the same as the
!interpolated velocity
!
!$OMP WORKSHARE
  Vxp(:,1)=Uxp(:,1)
  Vyp(:,1)=Uyp(:,1)
!$OMP END WORKSHARE
!
!For particles with inertia, we calculate the acceleartion
!
  do ist=2,nst
!$OMP WORKSHARE
    Axp(:,ist)=one_over_tau(ist-1)*(Uxp(:,ist)-Vxp(:,ist))
    Ayp(:,ist)=one_over_tau(ist-1)*(Uyp(:,ist)-Vyp(:,ist))
!$OMP END WORKSHARE
  enddo
!
!
endsubroutine get_particle_rhs
!!
!!---------------------------------------------------------
!!
subroutine get_velocity_at_particle_position(ux,uy,n1,n2)
  implicit none
  integer :: n1,n2
  double precision,dimension(n1,n2) :: ux,uy
  integer::i1,intx,inty,st
  real*8::xxp,yyp,xi,yj,xip1,yjp1,uuxp,uuyp

!!
  do st = 1,nst
!$OMP PARALLEL DO PRIVATE(i1,xxp,yyp,uuxp,uuyp)  SHARED(Xp,Yp,ux,uy,Uxp,Uyp)
    do i1 = 1,Nparticle
      xxp = Xp(i1,st); yyp = Yp(i1,st);
!
      call linear_interp_to_offgrid(xxp,yyp,n1,n2,deltax,deltay,ux,uuxp)
      Uxp(i1,st) = uuxp
      call linear_interp_to_offgrid(xxp,yyp,n1,n2,deltax,deltay,uy,uuyp)
      Uyp(i1,st) = uuyp
!
    enddo
!$OMP END PARALLEL DO
  enddo
!!
end subroutine get_velocity_at_particle_position
!!
!!------------------------------------------------------------
subroutine lambda_at_particle_position(GUxx,GUxy,GUyx,n1,n2)
  implicit none
  integer :: n1,n2
  double precision,dimension(n1,n2) :: GUxx,GUxy,GUyx
  integer:: i1,ist
  double precision :: xxp,yyp,gradup
!!
  do ist = 1,nst
!$OMP PARALLEL DO PRIVATE(i1,xxp,yyp,gradup)  SHARED(Xp,Yp,GUxx,GUxy,GUyx,dxuxp,dxuyp,dyuxp,Lamp)
    do i1 = 1,Nparticle
      xxp=Xp(i1,ist);yyp=Yp(i1,ist)
!
      call linear_interp_to_offgrid(xxp,yyp,n1,n2,deltax,deltay,GUxx,gradup)
      dxuxp(i1,ist)=gradup
      call linear_interp_to_offgrid(xxp,yyp,n1,n2,deltax,deltay,GUxy,gradup)
      dxuyp(i1,ist)=gradup
      call linear_interp_to_offgrid(xxp,yyp,n1,n2,deltax,deltay,GUyx,gradup)
      dyuxp(i1,ist)=gradup
      Lamp(i1,ist)=-dxuxp(i1,ist)*dxuxp(i1,ist)-dxuyp(i1,ist)*dyuxp(i1,ist)
    enddo
!$OMP END PARALLEL DO
  enddo
!
end subroutine lambda_at_particle_position
!!
subroutine calc_particle_diagnostic()
  write(*,*)Xp(1,1),Vxp(1,1),Uxp(1,1)
endsubroutine
subroutine get_particle_coordinate(pindex,stindex,xxp1,yyp1,vvxp1,vvyp1)
  integer :: pindex,stindex
  double precision :: xxp1,yyp1,vvxp1,vvyp1
  if((pindex.ge.1).and.(pindex.le.Nparticle).and.(stindex.ge.1).and.(stindex.le.nst)) then
    xxp1=Xp(pindex,stindex)
    yyp1=Yp(pindex,stindex)
    vvxp1=Vxp(pindex,stindex)
    vvyp1=Vyp(pindex,stindex)
  else
    write(*,*)'out of index'
    write(*,*)'Nparticle,nst=',Nparticle,nst
    write(*,*) 'pindex,stindex=',pindex,stindex
  endif
endsubroutine
!!-----------------------------
subroutine read_input_particle
  implicit none
  integer::i1,st,fno
  character*500 :: cnpar,formp
  character*80 :: fname1,fname2
!! Read the data for particle from input file
  write(cnpar,'(g8.0)') Nparticle
  formp = '('//trim(adjustl(cnpar))//'es16.3E2)'
  do st = 1,nst
    write(fname1,'(g8.0)')  st
    fno = 1000*st + 30
    open(unit=fno,file='input/tau'//trim(adjustl(fname1))//'/xptrack.out',status='old')
    read(fno,formp) (Xp(i1,st),i1=1,Nparticle)
    close(fno)
    fno = fno + 1
    open(unit=fno,file='input/tau'//trim(adjustl(fname1))//'/yptrack.out',status='old')
    read(fno,formp) (Yp(i1,st),i1=1,Nparticle)     
    close(fno)
  enddo
end subroutine read_input_particle
!%-----------------------------------
subroutine open_output_particle
  character*80 :: fname1,fname2
  integer::i1,fno,st
!!
  do st = 1,nst
    write(fname1,'(g8.0)')  st
    fno = 1000*st + 30
    open(unit=fno,file='tau'//trim(adjustl(fname1))//'/xptrack.out',status='unknown')
    fno = fno + 1
    open(unit=fno,file='tau'//trim(adjustl(fname1))//'/yptrack.out',status='unknown')
    fno = fno + 1
    open(unit=fno,file='tau'//trim(adjustl(fname1))//'/lamtrack.out',status='unknown')
    fno = fno + 1
    open(unit=fno,file='tau'//trim(adjustl(fname1))//'/dxux.out',status='unknown')
    fno = fno + 1
    open(unit=fno,file='tau'//trim(adjustl(fname1))//'/dxuy.out',status='unknown')
    fno = fno + 1
    open(unit=fno,file='tau'//trim(adjustl(fname1))//'/dyux.out',status='unknown')
    fno = fno + 1
    open(unit=fno,file='tau'//trim(adjustl(fname1))//'/vxint.out',status='unknown')
    fno = fno + 1
    open(unit=fno,file='tau'//trim(adjustl(fname1))//'/vyint.out',status='unknown')
    fno = fno + 1
    open(unit=fno,file='tau'//trim(adjustl(fname1))//'/uxint.out',status='unknown')
    fno = fno + 1
    open(unit=fno,file='tau'//trim(adjustl(fname1))//'/uyint.out',status='unknown')
  enddo
!!
end subroutine open_output_particle
!!
!! Write the data file for particle trajectory
subroutine write_output_particle
  implicit none
  integer::i1,fno,st
  real*8::x,y,xi,yj,xip1,yjp1
  character*500 :: cnpar,formp
!! 
  write(cnpar,'(g8.0)') Nparticle
  formp = '('//trim(adjustl(cnpar))//'es16.3E2)' 
!!
  do st = 1,nst
    fno = 1000*st + 30
    write(fno,formp) (Xp(i1,st),i1=1,Nparticle);	fno = fno + 1
    write(fno,formp) (Yp(i1,st),i1=1,Nparticle);	fno = fno + 1
    write(fno,formp) (Lamp(i1,st),i1=1,Nparticle);	fno = fno + 1
    write(fno,formp) (dxuxp(i1,st),i1=1,Nparticle);	fno = fno + 1
    write(fno,formp) (dxuyp(i1,st),i1=1,Nparticle);	fno = fno + 1
    write(fno,formp) (dyuxp(i1,st),i1=1,Nparticle);	fno = fno + 1
    write(fno,formp) (Vxp(i1,st),i1=1,Nparticle);	fno = fno + 1
    write(fno,formp) (Vyp(i1,st),i1=1,Nparticle);	fno = fno + 1
    write(fno,formp) (Uxp(i1,st),i1=1,Nparticle);	fno = fno + 1
    write(fno,formp) (Uyp(i1,st),i1=1,Nparticle);	fno = fno + 1
  enddo
    !!
end subroutine write_output_particle
!!
! Write the data file for particle lambda 
subroutine write_lambda_particle
  implicit none
  integer::i1,fno,st
  real*8::x,y,xi,yj,xip1,yjp1
  character*500 :: cnpar,formp
!! 
  write(cnpar,'(g8.0)') Nparticle
  formp = '('//trim(adjustl(cnpar))//'es16.3E2)' 
!!
  do st = 1,nst
    fno = 1000*st + 30
!    write(fno,formp) (Xp(i1,st),i1=1,Nparticle);       fno = fno + 1
!    write(fno,formp) (Yp(i1,st),i1=1,Nparticle);       fno = fno + 1
    fno = fno + 2
    write(fno,formp) (Lamp(i1,st),i1=1,Nparticle);      fno = fno + 1
!    write(fno,formp) (dxuxp(i1,st),i1=1,Nparticle);    fno = fno + 1
!    write(fno,formp) (dxuyp(i1,st),i1=1,Nparticle);    fno = fno + 1
!    write(fno,formp) (dyuxp(i1,st),i1=1,Nparticle);    fno = fno + 1
    fno = fno + 3
    write(fno,formp) (Vxp(i1,st)**2+Vyp(i1,st)**2,i1=1,Nparticle);      fno = fno + 1
!    write(fno,formp) (Vyp(i1,st),i1=1,Nparticle);      fno = fno + 1
!    write(fno,formp) (Uxp(i1,st),i1=1,Nparticle);      fno = fno + 1
!    write(fno,formp) (Uyp(i1,st),i1=1,Nparticle);      fno = fno + 1
  enddo
    !!
end subroutine write_lambda_particle
!!
subroutine write_trajectory_particle
  implicit none
  integer::i1,fno,st
  real*8::x,y,xi,yj,xip1,yjp1
  character*500 :: cnpar,formp
  character*80 :: fname1,fname2
!! 
  write(cnpar,'(g8.0)') Nparticle
  formp = '('//trim(adjustl(cnpar))//'es16.3E2)'
!!
  do st = 1,nst
    write(fname1,'(g8.0)')  st
    fno = 5000*st + 30
    open(unit=fno,file='tau'//trim(adjustl(fname1))//'/final_xptrack.out',status='unknown')
    write(fno,formp) (Xp(i1,st),i1=1,Nparticle);
    close(fno);
    fno = fno + 1
    open(unit=fno,file='tau'//trim(adjustl(fname1))//'/final_yptrack.out',status='unknown')
    write(fno,formp) (Yp(i1,st),i1=1,Nparticle);
    close(fno);
  enddo
    !!
end subroutine write_trajectory_particle
!!
!! Close the data file for particle trajectory
subroutine close_output_particle
implicit none
integer::i1,fno,st
!!
do st = 1,nst
  fno = 1000*st + 30
  do i1 = 1,10
    close(fno);
    fno = fno + 1
  enddo
enddo
!!
end subroutine close_output_particle
!! --------------------------------------
subroutine evolve_particles(dt)
  double precision :: dt
  call evolve_particles_euler(dt)
endsubroutine evolve_particles
!!-----------------------------
end module mod_part_interp
